# title:   Incentives for Non-participation: Absence in the United Kingdom House of Commons, 1997-2015
# content: script to reproduce results reported in the main paper, part 2
#          Figures 4-5 ("within MP" analysis)
# author:  Zoltan Fazekas
# note:    models take a lot of time to run
#          session info included in a separate document
#          assumes data file to be located in the same folder as script

# for RStudio warnings
# !diagnostics off

library("lme4")
library("tidyverse")
library("ggrepel")
library("merTools")
load("absence_within.Rdata")

# government to opposition model
gto <- glmer(absence ~ 1 + 
                elect_full + elect_second + in_second + 
                gvt_bill +
                factor(week_day_no) +
                const_maj_mm +
                dist_west_mm + sen_elapsed_mm +
                (1 + elect_full + elect_second + in_second | mp_name),
              data    = gvt_to_opp,
              family  = binomial(link = "logit"),
              control = glmerControl(optimizer = "bobyqa",
                                     optCtrl   = list(maxfun = 2e6)))
summary(gto)

#  opposition to government model
otg <- glmer(absence ~ 1 + 
                elect_full + elect_second + in_second + 
                gvt_bill +
                factor(week_day_no) +
                const_maj_mm +
                dist_west_mm + sen_elapsed_mm +
                (1 +  elect_full + elect_second + in_second | mp_name),
              data    = opp_to_gvt,
              family  = binomial(link = "logit"),
              control = glmerControl(optimizer = "bobyqa",
                                     optCtrl   = list(maxfun = 2e6)))
summary(otg)

# decompose for opposition party
otg_party <- glmer(absence ~ 1 + 
                elect_full + elect_second + in_second + 
                gvt_bill +
                factor(week_day_no) +
                const_maj_mm +
                dist_west_mm + sen_elapsed_mm +
                party +
                elect_full*party + 
                elect_second*party + 
                in_second*party + 
                (1 + elect_full + elect_second + in_second | mp_name),
              data    = opp_to_gvt,
              family  = binomial(link = "logit"),
              control = glmerControl(optimizer = "bobyqa",
                                     optCtrl   = list(maxfun = 2e6)))
summary(otg_party)

derivs1 <- otg_party@optinfo$derivs
sc_grad1 <- with(derivs1,
                 solve(Hessian, gradient))
max(abs(sc_grad1))
max(pmin(abs(sc_grad1),abs(derivs1$gradient)))

# create Figure 4
go_base <- gvt_to_opp %>% group_by(elect_full) %>% 
  summarize(elect_second = unique(elect_second),
            in_second    = unique(in_second)) %>% 
  ungroup() %>% 
  mutate(gvt_bill     = 1,
         week_day_no  = 0,
         const_maj_mm = mean(gvt_to_opp$const_maj_mm),
         dist_west_mm = mean(gvt_to_opp$dist_west_mm),
         sen_elapsed_mm = mean(gvt_to_opp$sen_elapsed_mm), 
         mp_name = "NA")
og_base <- opp_to_gvt %>% group_by(elect_full) %>% 
  summarize(elect_second = unique(elect_second),
            in_second    = unique(in_second)) %>% 
  ungroup() %>% 
  mutate(gvt_bill     = 1,
         week_day_no  = 0,
         const_maj_mm = mean(opp_to_gvt$const_maj_mm),
         dist_west_mm = mean(opp_to_gvt$dist_west_mm),
         sen_elapsed_mm = mean(opp_to_gvt$sen_elapsed_mm), 
         mp_name = "NA")

preds1 <- predictInterval(gto, data.frame(go_base), include.resid.var = FALSE,
                          level = 0.95, type = "probability")
preds1 <- cbind(go_base, preds1)
preds1$what <- "a) From government to opposition"

preds2 <- predictInterval(otg, data.frame(og_base), include.resid.var = FALSE,
                          level = 0.95, type = "probability")
preds2 <- cbind(og_base, preds2)
preds2$what <- "b) From opposition to government"

preds <- bind_rows(preds1, preds2)

ptime <- ggplot(preds, aes(x = elect_full, y = fit,
                           ymin = lwr, ymax = upr)) +
  geom_ribbon(alpha = 0.25) +
  geom_line(size = 1.1) +  facet_wrap(~what) +
  scale_y_continuous("", limits = c(0, 0.4)) +
  scale_x_continuous("", breaks = c(0, 0.5, 1),
                     labels = c("Beginning of\n 2005-2010",
                                "Elections", "End of\n2010-2015")) +
  theme_minimal() +
  scale_colour_manual("", values = c("grey50", "black"), 
                      labels = c("Opposition", "Government"),
                      guide = FALSE) +
  theme(text = element_text(family = "serif"), 
        panel.background = element_rect(fill = "grey98", colour = NA),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10),
        axis.title.y = element_text(size = 10, angle = 0, hjust = 1),
        axis.title.x = element_text(size = 10),
        plot.title = element_text(size = 10),
        panel.spacing = unit(3, "lines"),
        strip.text.x = element_text(face = "bold", hjust = 0),
        plot.margin = unit(c(1,1,1,1), "cm"))
ptime

# create Figure 5
og_base1 <- opp_to_gvt %>% group_by(elect_full) %>% 
  summarize(elect_second = unique(elect_second),
            in_second    = unique(in_second)) %>% 
  ungroup() %>% 
  mutate(gvt_bill     = 1,
         week_day_no  = 0,
         party = "Con",
         const_maj_mm = mean(opp_to_gvt$const_maj_mm),
         dist_west_mm = mean(opp_to_gvt$dist_west_mm),
         sen_elapsed_mm = mean(opp_to_gvt$sen_elapsed_mm), 
         mp_name = "NA")
og_base2 <- opp_to_gvt %>% group_by(elect_full) %>% 
  summarize(elect_second = unique(elect_second),
            in_second    = unique(in_second)) %>% 
  ungroup() %>% 
  mutate(gvt_bill     = 1,
         week_day_no  = 0,
         party = "LDem",
         const_maj_mm = mean(opp_to_gvt$const_maj_mm),
         dist_west_mm = mean(opp_to_gvt$dist_west_mm),
         sen_elapsed_mm = mean(opp_to_gvt$sen_elapsed_mm), 
         mp_name = "NA")

og_base <- bind_rows(og_base1, og_base2)
preds2 <- predictInterval(otg_party, data.frame(og_base), include.resid.var = FALSE,
                          level = 0.95, type = "probability")
preds2 <- cbind(og_base, preds2)
preds2$what <- "From opposition to government: party split"


p_party <- ggplot(preds2, aes(x = elect_full, y = fit,
                              ymin = lwr, ymax = upr,
                              linetype = party)) +
  geom_ribbon(alpha = 0.25) +
  geom_line(size = 1.1) +  
  facet_wrap(~what) +
  scale_y_continuous("", limits = c(0, 0.4)) +
  scale_x_continuous("", breaks = c(0, 0.5, 1),
                     labels = c("Beginning of\n 2005-2010",
                                "Elections", "End of\n2010-2015")) +
  theme_minimal() +
  scale_linetype_discrete("",
                          labels = c("Conservatives", "Liberal Democrats")) +
  theme(text = element_text(family = "serif"), 
        panel.background = element_rect(fill = "grey98", colour = NA),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10),
        axis.title.y = element_text(size = 10, angle = 0, hjust = 1),
        axis.title.x = element_text(size = 10),
        plot.title = element_text(size = 10),
        panel.spacing = unit(3, "lines"),
        strip.text.x = element_text(face = "bold", hjust = 0),
        plot.margin = unit(c(1,1,1,1), "cm"), legend.position = c(0.85, 0.2))
p_party